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REACTIVE BOUNDARY CONDITIONS AS LIMITS OF 
INTERACTION POTENTIALS FOR BROWNIAN 
AND LANGEVIN DYNAMICS* 

S. JONATHAN CHAPMANt, RADEK ERBAN*, AND SAMUEL A. ISAACSON? 

Abstract. A popular approach to modeling bimolecular reactions between diffusing molecules 
is through the use of reactive boundary conditions. One common model is the Smoluchowski partial 
adsorption condition, which uses a Robin boundary condition in the separation coordinate between 
two possible reactants. This boundary condition can be interpreted as an idealization of a reactive 
interaction potential model, in which a potential barrier must be surmounted before reactions can 
occur. In this work we show how the reactive boundary condition arises as the limit of an interac¬ 
tion potential encoding a steep barrier within a shrinking region in the particle separation, where 
molecules react instantly upon reaching the peak of the barrier. The limiting boundary condition 
is derived by the method of matched asymptotic expansions, and shown to depend critically on the 
relative rate of increase of the barrier height as the width of the potential is decreased. Limiting 
boundary conditions for the same interaction potential in both the overdamped Fokker-Planck equa¬ 
tion (Brownian dynamics), and the Kramers equation (Langevin dynamics) are investigated. It is 
shown that different scalings are required in the two models to recover reactive boundary conditions 
that are consistent in the high friction limit (where the Kramers equation solution converges to the 
solution of the Fokker-Planck equation). 
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1. Introduction. Let X{t) G il C and V(t) G denote the stochastic 
processes for position and velocity at time t of a molecule moving in the d-dimensional 
domain (where d gN) according to the Langevin dynamics (LD): 

dX{t) = V(t) dt, 

,_ ( 1 . 1 ) 

dV{t) = -{(3V{t) + D/3W(p{X{t)))dt + pV^dW{t), 

where /3 is the friction constant, D is the diffusion constant, tp : 12 —)■ M is the potential 
(to be specified later) and W(t) denotes a d-dimensional Brownian motion. Here /3 is 
assumed to have units of ‘per time’, D units of ‘distance squared per time’, and tp is 
assumed to be non-dimensional. In physical units, the potential energy of a molecule 
at X is then k-BTtf{x), where fee is the Boltzmann constant and T is the absolute 
temperature. If m denotes the mass of the molecule, the Einstein relation gives that 
mD[i = k-^T, so that excepting the noise term, (1.1) follows from Newton’s second 
law of motion. 

Passing to the overdamped limit /? —> cx) in (1.1), we obtain the Brownian dy- 
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namics (BD) model for X{t) as follows: 


dX{t) = DVip{X{t))dt + V^dW{t). 


( 1 . 2 ) 


In the special case that ip = 0 the molecule simply moves by Brownian motion, 


dX{t) = VwdW(t). 


(1.3) 


Equation (1.3) is a popular description for the movement of molecules within cells. 
It has been used to model spatial transport in a number of computational pack¬ 
ages for simulating intracellular processes, including Smoldyn [2,3,30], GFRD [34,35] 
and FPKMC [27,28] . A common approach for then coupling molecular interactions 
(diffusion-limited reactions) to (1.3) in these packages is to postulate that reactions 
can occur by one of several possible mechanisms if the corresponding reactants are 
sufficiently close [1,16[. 

The LD model (1.1) with tp = 0 provides a more microscopic description of diffu¬ 
sion than the BD model (1.3). It computes both position and velocity of a molecule by 
assuming that the molecule is subject to a normally distributed random force during 
each time increment. In particular, LD can be considered as an intermediate descrip¬ 
tion between detailed molecular dynamics (MD) simulations and BD simulators [12] . 
Typical full-atom MD simulations use time steps of the order of 1 fs = 10“^® s [24], 
while Smoldyn discretizes the equation (1.3) with times steps ranging from nanosec¬ 
onds to milliseconds, depending on a particular application [30] . Stochastic descrip¬ 
tions which compute both position and velocity of diffusing particles, including LD, 
are applicable on intermediate time scales [12,13]. 

One advantage of Smoldyn or similar BD packages is that they can simulate whole¬ 
cell dynamics. BD models based on equation (1.3) have been applied to a number 
of biological systems including signal transduction in E. coli [26], actin dynamics 
in filopodia [18], the MAPK pathway [34] and intracellular calcium dynamics [10[. 
In these applications, the positions of all diffusing molecules are updated according 
to (1.3) and the distances between each pair of possible reactants (for bimolecular 
reactions) are calculated. Each reaction then occurs (with a given probability) when 
the computed distance is smaller than a specified reaction radius (as in Smoldyn), or 
alternatively occurs when the computed distance exactly equals a specified reaction- 
radius (as in GFRD and FPKMG). To our knowledge, there is no established spatio- 
temporal simulator of intracellular processes based on the LD model. In order to 
develop one, one has to first investigate how bimolecular reactions might be described 
in the LD context. 

One possible way to implement bimolecular reactions in LD is to adopt the same 
approach as in BD. That is, the positions and velocities of a diffusing molecule would 
evolve according to the LD model (1.1) (with tp = 0) and each bimolecular reac¬ 
tion would occur (with a given probability) if the distance between two reactants is 
smaller than the reaction radius. However, as normally formulated this description of 
bimolecular interactions would not make use of a molecule’s velocity (as is available in 
LD). That is, the LD bimolecular reaction model would not provide any more physical 
detail than a BD model. 

In this work we step back from the normal BD bimolecular reaction model to 
the more microscopic reaction mechanism of a molecular interaction potential. The 
general potential forms we consider represent an irreversible bimolecular reaction as 
a molecule surmounting a steep potential barrier in the separation coordinate from a 
stationary target molecule, after which it enters an infinitely deep well (the “bound” 
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state). We show that the popular Smoluchowski partial-adsorption BD reaction 
model [22] can be derived in the simultaneous limit that the width of the barrier ap¬ 
proaches zero and the height of the barrier becomes infinite. Using the same potential 
model, we then examine what limiting reactive mechanism arises in the corresponding 
LD model, obtaining a specular reflection boundary condition. We conclude by show¬ 
ing that in the high friction limit, /? —> oo, the LD model with the specular reflection 
bimolecular reaction model converges back to the Smoluchowski partial-adsorption 
BD reaction model, consistent with the kinetic boundary layer studies of [7,8,23]. 
Our results demonstrate how an interaction potential model for a bimolecular reac¬ 
tion can be parametrized in either LD or BD models to be consistent with a BD model 
based on a partial-adsorption reaction mechanism. 

2. Problem Setup. We consider the movement of a molecule that can undergo 
an irreversible bimolecular reaction with a second stationary molecule, hereafter called 
the reactive target, which is modelled as a sphere of radius r^,. Let Br{0) C denote 
the d-dimensional ball of radius r centered at the origin. Then the reactive target 
is given as We assume that the diffusing molecule moves within the d- 

dimensional domain U which satisfies 

(0) \ (0)^ C ft, for rh<Li<oo. (2.1) 

Equation (2.1) defines U as the domain which is the exterior to the reactive target 
(sphere 9B,.j^(0)) and which includes all points which have a distance less than Li — 
Tb > 0 from the reactive target. A simple example of U satisfying (2.1), for which we 
can And some explicit solutions in Section 3.1, is given as 

n = R‘^\Br^{0). ( 2 . 2 ) 

This is a standard ansatz for deriving formulae for the reaction radius in BD descrip¬ 
tions, e.g. in the Smoluchowski or Doi models [1,16,32]. In most cases, we will further 
assume that U is bounded with smooth boundary, i.e. equation (2.1) is altered by a 
bound from above 

(^Bli(O) \ Bri,(0)^ C U C 5 ^ 2 ( 0 ), for rb < Li < L 2 < 00 . (2.3) 

A spherically symmetric example of domain U satisfying condition (2.3) is given as 

U = Bi(0) \ Brj,(0), for rb < L < 00 . (2.4) 

If dimension d = 1, then by the “surface of the ball” we mean the origin, i.e. ^Br^^ (0) = 
{0} and rb = 0. Then equation (2.4) reduces to the finite interval U = [0,L], which 
we use in our numerical examples. 

We assume that the molecule is adsorbed instantly upon reaching the surface of 
the reactive target, so that 

trajectory X{t) is terminated, if X{t) G ^Br^^(0). (2.5) 

In order to reach the surface, the molecule has to overcome a potential barrier, which, 
denoting r = j a: j, is given by 



r - rb\ 

£ J’ 


(p(x) = ip{r) 


for Cb < r < rb -I- e, 
for r > r\^ + e, 


( 2 . 6 ) 
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where e > 0 and </? > 0 are positive constants, and ip : [0,1] —>■ [0,1] is a smooth 
function with a (unique) global maximum at ip{0) = 1, and ip{l) = 0. These imply 

ip{z) < V'(O), 0 < z < 1, 


with 


dip 


( 0 ) < 0 . 


We assume that e in (2.6) is chosen sufficiently small so that e <C (Ti — rt,) for ri, and 
Li given in assumption (2.1). 

Instead of studying the Langevin equation (1.1) we shall work with the corre¬ 
sponding equation for the probability density that X{t) = x and V{t) = v, denoted 
p{x,v,t), which satisfies the Kramers equation (also called the phase-space Fokker- 
Planck equation) 


+ V -VxP = ■ [vp ++ P DV^p] , for X G ^,v G R'^, (2.7) 

where (resp. Vt,) denotes the gradient in x (resp. v) variable. Considering the 
overdamped limit (1.2), the corresponding equation for the probability distribution 
p{x,t) is given as the Fokker-Planck equation 


— = ZlVa, • [VxP-f pVa,<p], for a; e fl. (2.8) 

In this paper, we show that equations (2.7) and (2.8) with fully adsorbing boundary 
condition (2.5) on are in suitable limits equivalent to a diffusion process with 

partially adsorbing (reactive, Robin) boundary condition on dBr^(0), namely to the 
problem 

dp 

— = 77 A^p, for X Gfl, (2.9) 

with 

DVxPix) ■ — = Kp(x), for X G dBri^{0), (2.10) 


where K is & suitable Robin boundary constant (reactivity of the boundary) and 
x/r\^ is a unit normal vector to sphere dBj.^{0) at point x. Considering spherically 
symmetric fl given by (2.2) or (2.4), and a spherically symmetric initial condition, we 
can rewrite the Fokker-Planck equation (2.8) as 


dp 

~dt 


D d 
1 



dp 

dr 



( 2 . 11 ) 


where p{r, t) is the radial distribution function. Then the limiting Robin boundary 
problem (2.9)-(2.10) is given as 


dt dr \ dr J ’ 


dv 

with D -^{r\P) = Kp{r\P). 


( 2 . 12 ) 


The Robin boundary constant K in equations (2.10) and (2.12) is (together with D) 
an experimentally determinable (macroscopic) parameter. In this way, we are able to 
parametrize both BD and LD models using experimental data. 
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The rest of this paper is organized as follows. In Section 3, we investigate the 
BD description (1.2). To get some insights into this problem, we first consider the 
specific case of a linear interaction potential and spherically symmetric domain (2.2) 
for d = 3. In this case, we can explicitly solve the corresponding Fokker-Planck 
equation as shown in Section 3.1, and prove the convergence of this solution, as e —> 0 
and ^ —>■ 00 , to the solution of a model involving a reactive Robin boundary condition. 
We then continue in Section 3.2 with an asymptotic analysis, in the same dual limit, 
of the full BD model (1.2) in general domain (2.3). In Section 4, we investigate the 
LD model (1.1) and derive a boundary condition in the limit e —>■ 0. This boundary 
condition is then used in Section 5 to connect the LD model to a BD model with 
Robin boundary condition in the dual limits of high friction, ^ ^ oo, and large 

potential barrier, ip —> oo. Numerical examples supporting our analysis are provided 
in Sections 4.1 and 5.1. We conclude with discussion in Section 6. 

3. Brownian Dynamics. In this section we consider the overdamped problem 
in which the particle moves by BD, i.e. its position X{t) evolves according to equa¬ 
tion (1.2) with the interaction potential given by (2.6). 

3.1. Simple three-dimensional example with explicit solution. Before 
investigating the general problem in as a warm-up we first consider a simplified 
special case to gain insight into how to recover the Robin boundary condition (2.12) 
from an interaction potential. We consider the spherically symmetric domain (2.2) 
for d = 3. The steady-state spherically-symmetric Fokker-Planck equation (2.11) is 
then given by 


J.2 Qj. 


dp dip 


for r\y < r < oo, (3.1) 


where the fully adsorbing boundary condition (2.5) implies a Dirichlet boundary con¬ 
dition at r = Tb, i.e. p{rh) = 0. We consider constant concentration p^o > 0 far away 
from the reactive surface, i.e.^ 


limp(r)=Poo- (3.2) 

r—>-oo 

We also assume that the potential p{r) is given by (2.6) where ijj : [0,1] —)■ [0,1] is 
the linear function 


Then 


ipiz) = 1 — z. 


/ X _ J - 

= 


'F 


rh < r < Tb + s, 
r\^+ e < r. 


(3.3) 


Substituting into (3.1), and making use of the boundary conditions at r = rb and 
r = oo, we can solve piecewise to obtain a solution on (rb,rb + £) and a solution on 


^Of course, p as defined here is not a probability distribution, since it does not integrate to unity 
on the infinite domain r > rb- The problem as stated arises from a limiting process in which p is 
re-scaled appropriately as the domain size tends to infinity. 
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(rb + e,oo). On each interval the solution is defined up to an unknown constant. By 
enforcing continuity of the flux at r = rb we can eliminate one constant to obtain 



'a r \ exp 

1 , 

1_ 

p{r) = < 

Jr, 

Poo 7 

k r 

L £ J 


ds, 


for Th < r < r\y + £, 
for r > rb + e, 


(3.4) 


where the last constant, A, is determined by the continuity of p(r, t) at r = rb + £. 
This gives 


(rb + £) Poo 

1 + {rh + £)/3(e,(^)’ 


with 


I3{£,‘p) 



exp[(^(l-s)] 

ds 

(rb +es)^ 


(3.5) 


To recover the Robin condition (2.12) as e —)■ 0 and ip —)■ oo, we need 


dp 

lim lim D — (r) — Kp(r) 

_e->-0 r->-(ri,+e)+ dr 

(p—^OO 


= 0 , 


which is equivalent to 


lim 
£ —^0 
)-oo 


D + K{rh+e) 

(rb + s) (l + (rb + £)/3{£, p)) 


= K, 


or simplifying, 


lim (3{£, (p) = 


D 


(p—>-oo 


We note that 


0<l3{£,p)<^ / exp[<p(l - s) ] ds = 


e(exp[(/3] - 1) 


'b -'o 


rip 


(3.6) 


(3.7) 


Therefore, f3{£,p) will have a finite limit if ee'x.p[p]/p has a finite limit. This moti¬ 
vates the choice 


£ = 


D p exp[ —p] 
K ' 


Using (3.5), we have 


I3{£,p)-^ [ exp[</5(l - s)]ds = e[ exp[(^(l - s) ] ( ^-^ 

rt Jo Jo \ (rh + 


1 

es)2 rl 


(3.8) 


ds 


< 


De (2rb -I- e) 


Krl 


(3.9) 
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which converges to zero as e —>■ 0. We therefore conclude that 

- e /■! _ 

lim B(e,ip) = lim / exp[ (p(l — s) 1 ds 

£->■0, v-^-oo, E-fO, ’’b-to 

where e and Tp are where e and Ip are 

related by (3.8) related by (3.8) 

D^exp[-^] D 

= hm -—^- / exp (/?(!-sj ds = 

Kr^ Jq Kr^ 

so that we get (3.6). In particular, we recover the Robin boundary condition with 
Robin constant K. 

The steady-state solution to the limiting Robin boundary condition problem (2.12) 
with (3.2) is given by 


p{r) = po 


1 - 


Th 


1 -L 


D 

Krh 


for Th < r < 00 . 


We now examine the error between this limit and the solution (3.4) for two cases: 
r\y < r < + e and r > + e. For the latter we have 


|p(r) -p(r)| = 


We can estimate that the denominator is greater than r. Consequently, 


Poo 

e 

,1 + /3(e, V^)) 

+ rl 



r 

1 + (Lb -L£)/3(£,(^)) 




\p{r) -p{r)\ < 


Poo e 


1-f 


D 

Kry, 




Poo rl 




Using (3.7) and the definition of e (equation (3.8)), we have f3{e,'ip) ~ 0(1) as e —> 0. 
Using (3.9), we conclude that the second term on the right hand side is 0(e/r) as 
e —>■ 0. Thus, we conclude 


\p{r) — p{r)\ = O = 0{e), for r>ri^+e. 

On the other hand, we have 

Pqq D 

sup \p{r) - p{r) I > |p(rb) - p(rb) | = p(rb) = „ „ 

rG(rb,Tb+£) K.T}q U 


(3.10) 


> 0 . 


The maximum error between the two models is therefore 0(1) as £ —0 (where Tp 
satisfies (3.8)), illustrating the non-uniformity of the error within the region where 
the potential interaction is non-zero. 

Note, while the maximum norm (i.e. norm) of the difference between solutions 
does not converge on (rb,rb -I- £), both p{r) and p(r) are uniformly bounded on 
(Lb, Lb + e) since 

sup p{r) < A(3{e,p) < 

'uG(rb,rb+£) 

sup p{r) < Poo- 
i'G(r-b,rb+e) 
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As such, using that the maximum norm error is 0{e/r) on (rb + e,oo) by (3.10), we 
find that the g-norm converges for any 3 < q < oo, 

\\p{r) - p{r)\\q = \p{r) - p{r)\‘‘ dr^ = O (e?) . 

Here, the lower bound g > 3 is an artifact of our working on an unbounded domain, 
which requires integrability of |p(r) —p(r)|'^r^ on (rb,oo). If our domain was bounded, 
i.e. given by equation (2.4), and we used the Dirichlet boundary condition that 
p{L) = Poo, we would expect a similar estimate to hold for all 1 < g < oo. 

In summary, we have found that in the special case of a linear potential barrier, 
by taking the height of the barrier p —>■ oo and then choosing the width of the barrier, 
e, to decrease exponentially in the height (i.e. according to equation (3.8)), we can 
recover the solution to the diffusion equation with Robin boundary condition. We 
may therefore interpret the Robin boundary condition model for bimolecular reactions 
between two molecules as approximating an underlying interaction potential, in which 
the two molecules must surmount a steep potential barrier before entering a bound 
state represented by a deep well. 

Remark. In [31] it was shown how the Robin constant can be chosen to give the 
same diffusion-limited reaction rate in two steady-state spherically-symmetric models 
both including the same fixed interaction potential, with one using a zero Dirichlet 
boundary condition at rb (as in (3.1)), and the other using a Robin condition at rb+e 
for any fixed e > 0. The argument in [31] can be modified to match diffusion-limited 
reaction rates in the two models considered in this section (i.e. (3.1) with p(rb) = 0, 
and the steady-state diffusion equation with Robin boundary condition satisfied by 
p(r)). It is satisfying to note that choosing K to match the diffusion limited rates 
in this manner also gives the relationship (3.8) between K, e and p. We therefore 
conclude that choosing K to match the diffusion-limited reaction rates of the two 
models of this section is sufficient to give convergence of p{r) to p{r) as e —>■ 0 and 
p —t 00. 

3.2. Asymptotics for a general BD model with a general interaction 
potential. To analyze equation (2.8) for arbitrary d G N with general potential p(r) 
in the form (2.6), we follow the approach of Pego [29] and Li [25] . We use the method 
of matched asymptotic expansions in a general bounded domain n which satisfies 
(2.3) and has a smooth boundary. First, we consider an inner solution in a boundary 
layer near dBr^,(0). Let 


z = 



(3.11) 


denote the stretched distance of a point x from the reactive boundary. We define an 
inner solution p{z, x, t). The second argument should actually he x — x/ \x\, but we 
follow the method of [29] and instead assume that 


p{z,Xx,t) =p{z,x,t), 


for all A > 0 (i.e. the length of x is accounted for by z). We note the identities that 


X/ X ' ^ — 


d-1 

r 



d-1 
rb + ez ’ 


(3.12) 
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In the new coordinate system 


^ ^ , I dp ^ 

s oz 

A A~ I d — 1 dp 

A^p A^p H-;-— 

e rh + ez oz 


1 d'^p 
dz'^ ’ 


where we have used (3.12) and that 


- ^ dp . , 


= 0 , 


y=x 


as p is assumed constant when the second argument is varied along the x direction. 
From (2.8) we find the inner problem (for 0 < z < 1) 


dp 

dt 


= D 


AxP- 


1 d — 1 dp 
e rh + ez dz 


1 d'^p 

£.2 Q^2 


ip dp dll) 
dz dz 


ip d — 1 d'i/) 

P 


+ ^P 


dV 


e Tb + ez dz dz^ 


We now expand the inner and outer solutions in e and denote the leading-order terms 
by pQ and pq respectively. The leading-order behavior of the inner solution then 
satisfies 

d'^Po dpodtp ^ ^ d^ 

Note that we will find below that in the distinguished limit p depends on e in a 
logarithmic manner (see equation (3.13)) so that retaining p here is consistent with 
neglecting terms of 0(e). Solving this equation and using the Dirichlet boundary 
condition at z = 0, we find 

Po{z,x,t) = A{x,t) / exp[p(V'(z') - V'(^))] dz', 

Jo 

where A(x,t) is an unknown constant, also satisfying the dilation property 
A{XxA) = A{x,t), for A > 0. 

The leading order term of the outer solution, po, satisfies the diffusion equation (2.9) 
away from the reactive boundary. Our matching conditions are 


x^dB 


lim pAxA) = AmpQ{z,x,t) = A{x,t) / expr(p^(z')l dz', 

asr-jo) Jo 


lim Wpo{x, t) ■ X = 


lim 


>9^0 1 

= iim - 


X^dBr^ (0) 


x^dBrAO) dr z->-l £ 


dpo , dp) 
oz dz 


A(x, t) 


Combining these two equations, for x g dBr^{0) we find 

i ^ r^^PO, DPoix,t) Dp 

DVxV\^) ' — = ^- 

dr £ exp[jpp>{z')] dz' £exp(<p) 


dtp 

dz 


( 0 ) 


Poix,t), 


where the last approximation is obtained using Laplace’s method in the limit p —>■ oo. 
Thus we recover the desired Robin boundary condition (2.10) if p —>■ oo and £ —>■ 0 
such that 


D ip 

K exp(^) 


dtp 

Jz 


( 0 ) 


e = 


(3.13) 
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In summary, we have used the method of matched asymptotic expansions to ex¬ 
amine scaling limits for a general set of potential interactions between two particles. 
The potentials were assumed to be short-range, forming a high barrier as the separa¬ 
tion between the molecules approaches a fixed “reaction-radius”, and then a deep well 
once this barrier is surmounted. In the limit that the height of the barrier, ^ > oo, 
and the width of the barrier, e, decreases to zero exponentially in the height, we re¬ 
cover the solution to the diffusion equation with Robin boundary condition. We may 
therefore interpret bimolecular reactions between two molecules modeled with a Robin 
boundary condition as an approximation to one of many possible underlying poten¬ 
tial interactions. These interactions are characterized by the two molecules needing 
to surmount a steep potential barrier before entering a bound state represented by a 
deep well. 

Remark. If e is chosen to approach zero slower than (3.13), then we recover a zero 
Neumann boundary condition at the reactive surface, 

= 0, for X S (3-14) 

Likewise, if e is chosen to approach zero faster than (3.13), then we recover the zero 
Dirichlet boundary condition that 

p{x,t)= 0 , for a; G 9Rrb(0)- 

4. Langevin Dynamics. We now consider the LD model (1.1) in a bounded 
d-dimensional domain satisfying (2.3). The reactive target is again taken to be the 
surface of the d-dimensional sphere of radius rb about the origin. It is again assumed 
that the molecule is adsorbed instantly upon reaching the surface of the sphere, i.e. we 
consider the boundary condition (2.5) together with spherically symmetric interaction 
potential p given by (2.6). We leave unspecified any boundary condition on 912 \ 
9i?rb(0)) it is not needed in the following analysis. 

Instead of studying the Langevin equation (1.1) we work with the corresponding 
Kramers equation (2.7). Since x/r-h is the normal to 9i?r.,,(0) at x, the adsorbing 
boundary condition (2.5) means that the Kramers equation (2.7) is coupled with the 
Dirichlet boundary condition 

p{x,v,t) =0 for a; G dBr^,{0) and v ■ x > 0, (4.1) 

and the unspecified boundary condition on dn\dBr^^(0). We are interested in various 
limits of (2.7) as e —>■ 0, ^ (X) and /3 —>■ (X) in which the interaction potential can 
be approximated by an appropriate reactive boundary condition. In the BD models 
of Section 3 our goal was to derive the widely-used Robin boundary condition. To 
our knowledge, in the LD model we now investigate there is no standard reactive 
boundary condition for bimolecular reactions. Therefore, we wish to see what, if any, 
reactive boundary condition arises when considering similar limits of s and ip. 

As in Section 3.2, to study these limits we will match an inner solution within an 
0{e) boundary layer about dBr^^(0) to an outer solution when far from the reactive 
boundary. Using the same notation as in Section 3.2, we denote by z the stretched 
distance from the boundary, given by (3.11). We also introduce a re-scaled radial 
velocity, 

Vr V ■ X 
V = — = — 


V; 


V, 



REACTIVE BOUNDARY CONDITIONS AS LIMITS OF INTERACTION POTENTIALS 11 


where x = a:/ |a;| is a unit vector in the x direction, Vr = v ■ x is the radial velocity 
into n from and Vr is a scaling constant in the radial velocity that will be 

specified later. In the inner region we denote the density by p(z, x, v, v,t), where we 
assume that p is constant whenever x is varied in the radial direction, and when the 
component of the velocity, v, in the radial direction is varied. That is 

p{z, X, I/, V, t) = p{z, Xx, i/, V, t), for A > 0, 

p{z, X, I/, V, t) = p{z, x,h',v + ax, t), for a e K. 


In addition to the identities (3.12), we note that 

V vx 

'SI — 


Vr{rh + £z) Th + ez' 


(4.2) 


Using (3.12) and (4.2) we find the derivative operators in (2.7) transform in the new 
coordinates to 


^ ^ . I dp ^ 

Va;P+ -ZrX + 

e az 

^ ^ _ 1 dp ^ 

V^p+ —z^X, 

Vr dv 


V X 


dp 

Pr (Lb + £z) + ez J dv’ 




1 d^p 

Vr dv"^ ’ 


where we have used that 


dp 


X ■ X/^,^{z,x,v,v',t) 
ov 

Using (2.6), equation (2.7) can be transformed to 


= 0 . 


dp ^ VrV dp f V ■ V VrV^ \ dp 

dt ^ a:P+ ^ dz~^ \Vr{rh+ £z) r\y + £Z) dv 

o( i~ TT ~ P^P p D dip dp PO d‘^p\ 

jdldp + v ■ V^p + v^ + U— -T ^ ^ PO/X^pd- ■ 

\ ov VrS dz ov Vr OV‘^ / 


We consider an asymptotic expansion of the inner solution, p (in the boundary layer 
about dBr^,{0)) as e —>■ 0 with all other parameters held fixed, 

p ~ + p^^^e + pP^^e^ + _ 


Similarly we also consider an expansion of the outer solution, valid far from dBr^,{0), 
for which we will abuse notation and also denote it by p, 

p ~ +_ 


At leading order 0{e ^) we find that the inner solution satisfies 

5p(o) p P D dip dp^^’^ 

^ Ai-^— A-= 0- 

OZ Vr dz OV 


To simplify this equation, we let 

Vr = y^pP D. 


(4.3) 


( 4 . 4 ) 
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velocity at x = e 


Fig. 4.1: (a) Characteristic curves of (4.5) for the linear potential (3.3) are shown as 
blue lines. The curve z = i/^/2 is drawn as a thicker black line, dividing the upper 
half plane into a region where the moving particle will escape over the barrier (v < 0, 
with z < 12), a region where all trajectories are reflected (z > 12), and a region 

where p^^^ = 0 with all trajectories originating on the z = 0 axis (v > 0 and z < 12.) 

because of boundary condition (4.6). 

(b) Probability of escaping through the left boundary as a function of the incoming 
velocity for e = 0.03 (blue line), e = (green line) and e = 10“^ (red line). We 

compute the trajectories using (4.9)-(4.10) for At = 10“^, D = 1, L = 1, (j = 10^ and 
if = 2.593. The vertical dashed line at —\/2Tp (3 D separates the two cases of boundary 
condition (4.7). 


This choice emphasizes large velocities, for which we expect the particle to be most 
likely to escape over the (reactive) potential barrier. Substituting (4.4) into (4.3), we 
obtain 


dz dz dv 

The boundary condition (4.1) implies 

p^^\d,x,v,v,t) = 0, for 1 /> 0. 


(4.5) 


(4.6) 


The solution to (4.5) is constant along the characteristic curves 

-V'(^) = y + 

for C an arbitrary constant. These curves are shown in Figure 4.1(a) for the linear 
potential (3.3). The curve 


iP' 

1 - V'(^) = y 

divides the half-plane into three regions. In the region where ly > 0 and 1 — if^z) < 
(2, p'T'l is zero due to the boundary condition (4.6). This region corresponds to 
particles originating at the boundary and moving into the domain. Since particles 
are only adsorbed at the boundary, and not emitted, we find that the solution is 
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zero throughout the region. Where < 0 and 1 — < z/^/2 the solution will be 

determined by matching to the outer solution. This region corresponds to particles 
entering from outside the boundary layer with sufficient velocity to escape over the 
potential barrier, and thereby exit through the reactive surface at z = 0. Finally, 
in the region 1 — z/’(z) > z^^/2, trajectories with v < 0 are reflected in a symmetric 
manner. For points with z/ < 0 the solution will again be given by matching to the 
outer solution, while for those with zz > 0 the solution is determined by reflection. 
This region corresponds to particles that enter from outside the boundary layer with 
insufficient velocity to escape over the potential barrier. These particles are then 
reflected, moving back out of the boundary layer. Note, in the original variables the 
curve separating these three regions is given by 


1 — 'ijj 



{v ■ xY 

2tpJd ~ 2^Jd' 


For X G we match the inner solution as z approaches the edge of the bound¬ 

ary layer to the outer solution as the particle’s spatial position, y, approaches x. That 
is, we require 

lim X, V, V, t) = lim v, t). 

2—>-1 y^x 


Using that '0(1) = 0, this matching condition implies the outer solution satisfies the 
following reactive boundary condition for Vr > 0 and x G (0) 


p^°\xyV,t) 


(a:, V — 2vr x,t), 0 < Vr < 

0, Vr > y/2lpl3 D. 


(4.7) 


When the moving particle reaches the reactive boundary with radial velocity in the 
outward direction greater than ^/2ip /? D it leaves the domain (i.e. undergoes reac¬ 
tion). In contrast, when the particle reaches the boundary with a slower radial velocity 
in the outward direction it is reflected back into the domain along the direction of 
the normal at the point where it hit the reactive boundary. This “specular reflection” 
boundary condition is also obtained in the corresponding deterministic Newtonian me¬ 
chanics model. We demonstrate this explicitly for a simple one-dimensional example 
in Appendix A. 

A version of this boundary condition, given in terms of an arbitrary threshold 
velocity, is assumed in the kinetic boundary layer investigations of the one-dimensional 
Kramer’s equation in [7,8], and the 3D spherically-symmetric steady-state Kramer’s 
equation in [23]. It is also (briefly) mentioned in [23] that the specific threshold of 
2^/30 we derive is what one might impose across an interface where the potential is 
discontinuous with a jump of size Tp (again, for the 3D spherically-symmetric steady- 
state Kramer’s equation). Our asymptotic analysis shows hows this specific threshold 
velocity arises in the general Kramer’s equation as the limit of a shrinking potential 
boundary layer bordering a Dirichlet boundary condition. We obtain an effective jump 
in potential at the reactive boundary, as opposed to an interface within the domain. 
Contrast this to the limit of the BD problem from Section 3.2, where in taking e —>■ 0 
with Tp fixed the influence of the potential is completely lost (e.g. a zero Dirichlet 
boundary condition is recovered). 

We therefore conclude that the LD model with the interaction potential is in the 
limit that e —^ 0 equivalent to solving the (zero-potential) Kramers equation 

dj) 

-£ + v- WyrP = /? v„ • [vp + l3D'G/y,p ], for a; G e 


(4.8) 









14 


S. J. CHAPMAN, R. ERBAN, S. A. ISAACSON 


with the specular reflection reactive boundary condition (4.7) on (0) and what¬ 
ever boundary condition was imposed on (90. That is, in the limit that the width 
of the potential approaches zero, with the barrier height held flxed, we And that the 
potential can be approximated by a velocity threshold boundary condition. Here par¬ 
ticles moving sufficiently fast relative to the height of the potential barrier undergo 
bimolecular reactions when reaching the reactive boundary, while those moving too 
slow are reflected back into the domain. This result should be applicable for general 
short-range potential interactions that form a high barrier as the separation between 
two molecules approaches a flxed “reaction-radius”, and then a deep well once this 
barrier is surmounted. In contrast to the diffusive case, taking the barrier height 
Ip ^ oo (with (3 flxed) leads to a complete loss of reaction; all particles reaching the 
reactive boundary are simply reflected back into the domain. 

4.1. A numerical example showing recovery of boundary condition (4.7) 
as e —>■ 0. We consider the LD model (1.1) for d = 1 and the linear potential (3.3). In 
the one-dimensional case, our computational domain is interval H = [0, L\. We choose 
a small time step At and compute the position X{t + At) and the velocity V{t + At) 
from the position X (t) and the velocity V(t) by 

X(t + At) = X{t) + V(t)At, (4.9) 

Vit + At) = V{t) - (3 V{t) At + X[o,e] {X{t)) At + j3'j2DAt (4.10) 

where X[o,£] : R —> {0; 1} is the characteristic function of the interval [0,e] and ^ 
is a normally distributed random variable with zero mean and unit variance. We 
implement adsorbing boundaries at both ends a; = 0 and x = L oi the simulation 
domain [0, L] by terminating the computed trajectory whenever X{t) <0 or X{t) > L. 

In this paper, we are interested in understanding the dependence of the behavior 
of the LD model (1.1) on its parameters e, Ip and /3. In particular, we choose the 
values of other parameters equal to 1, namely 

D = L = 1. (4.11) 

In this section, we are interested in the limit e —>■ 0. We want to illustrate the 
boundary condition (4.7). Therefore we fix the values of Ip and /3 and simulate the 
LD model (1.1) for different values of e. 

For each value of e, we compute many trajectories according to (4.9)-(4.10), 
starting from the middle of the domain, i.e. A(0) = L/2, with initial velocity 1^(0) 
sampled from the normal distribution with zero mean and variance j3 D. Whenever a 
particle enters the region [0,e], we record its incoming velocity. Then we follow its 
trajectory in the region [0, e] and record one of two possible outputs: 

(1) the particle leaves fl through its left boundary (i.e. X{t) < 0); or 

(2) the particle returns back to the region (e, L] of domain il (i.e. X{t) > e). 

The fraction of particles which left the domain D through its left boundary as a 
function of the incoming velocity is plotted in Figure 4.1(b). This curve can be 
interpreted as the probability that the particle escapes over the potential barrier, 
given its incoming velocity. We also plot the vertical dashed line at (3 D in 

Figure 4.1(b). This threshold separates the two cases of boundary condition (4.7). 
We observe that the probability of escape converges to the step function as e —>■ 0, 
i.e. we have numerically confirmed boundary condition (4.7). 

We will return to this example in Section 5.1 when we study the convergence of the 
LD model to the diffusion process with a Robin boundary condition. In particular, we 
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use values /3 = 10^ and Tp = 2.593 in Figure 4.1(b). This choice of Ip will be explained 
in the following section and then used again in one of our simulations presented in 
Section 5.1. 

5. From Langevin Dynamics to Brownian Dynamics. We now study the 
overdamped limit of the LD model (1.1). We wish to show that in the overdamped 
limit where /3 —> cx), taking ^ oo in a /3-dependent manner will recover a Robin 
boundary condition for the limiting diffusion equation. To study the /3 —)■ cx) limit, 
we extend the asymptotic analysis of the one-dimensional Kramers equation in [14] 
to the d-dimensional (zero-potential) Kramers equation (4.8). The kinetic boundary 
layer studies in [7, 8, 23] previously investigated the relationship between the velocity 
threshold in the specular reflection boundary condition and effective adsorption rate, 
K, in the Robin condition. Our approach here differs from those studies, which 
primarily used numerical solutions of truncated moment equations or basis function 
expansions to estimate empirically determined formulas for the Robin constant, K. 
We focus on deriving an explicit formula relating how the potential barrier height, Tp, 
should be chosen as /3 —>■ c» to recover a specified Robin constant. 

We begin by re-scaling velocity as u = ^/]3 t] and let f{x, t], t) = p(x, v, t). Sub¬ 
stituting into (4.8), we obtain 

|[ + V^J7 -V./ = /3V^-]r,p + DV^p]. (5.1) 

We expand / in powers of /3“^/^ as 

f{x, T 7 , t) - fo{x, r],t) + ^ fi{x, ^ f 2 {x, r],t) + .... 

Substituting into (5.1), we find 


V,, • [nfo+D\/rjfo] = 0, 

(5.2) 

'^v ■ [nfi + A’V^/i] = T] ■ Vx/o, 

(5.3) 

• ]»7/2 + D V^/ 2 ] = J7 • V,/i + ^. 

(5.4) 


Implicit in these equations is the assumption that we are interested in timescales for 
which 


We therefore interpret (5.2) as implying that the velocity distribution component 
of /o relaxes to equilibrium on a faster timescale than t. Denote by r this faster 
timescale. As discussed in the introduction to ]9], up to a normalization constant 
there is a unique solution to (5.2) that corresponds to the equilibrium solution of the 
fast-timescale, time-dependent equation 


dfo 

dr 


= V^-]r7/o-bDV^/o]. 


This equilibrium solution is then 


fo{x,r],t) = g{x,t) exp 


1 ^ 1 ^ 

2D 
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where g{x,t) is independent of rj. We similarly find that the general solution to (5.3) 
is given by 


1^1 

2D 


fi{x,T],t) = (^-V^g{x,t) ■'n + Cix,t)j exp 
where ^{x,t) is also independent of r/. Substituting these into (5.4) we see that 

V^-[r,/2+i^V^/2]= ^ 


dxidxj 
i,j=i ■> 



1 ^1^ 

exp 

2D 


Integrating both sides of this equation for all r/ G and using that f 2 {x,ri,t) -G 0 
and \Wrjf 2 {x,ri,t)\ —>■ 0 as \ri\ —>■ oo, we conclude that g{x,t) satisfies the diffusion 
equation 

dg 


dt 


= Dl^^g. 


The probability density that the particle has position x at time t and has not reacted 
is given by 


u{x,t) = / p{x,v,t)dv, 

Jr<^ 


implying that to leading order u is proportional to g. As such, we expect as /3 —)■ cx) 
that u satisfies the diffusion equation. 

We now show that g, and hence u, satisfies a Robin boundary condition as /3 ^ oo 
when ip is chosen to approach infinity in a /3-dependent manner. To leading order in 
P, the outward flux through a point, x G is 

J{x,t): = — {v ■ x) p{x,v,t) dv 
JR‘‘ 

- [ {v-x) fo{x,v/^/p,t) + ^ fl{x,v/^,t) dt! 

jR-i L VP 

= {2ttD [p V^gix, t) ■x'^, 
where as in previous sections x = x/ \x\. Let 

TZ[x) = |r; 1 1 ! • £ < —\j2Tp(3 d'^ = |t! 


(5.5) 


Vr < — 




denote the set of velocities at which particles may escape through the reactive bound¬ 
ary. The specular reflection boundary condition (4.7) implies 


J{x,t) = — / [v ■ x)p{x,v,t) 

~ / (^ • *) 

J'R.ix) 


dt! 


in{i 
= (2 7rD/3)‘^/2 


fo{x,v/^,t) + fi{x,v/y/p,t) 


dt! 


exp[-v3] g{x,t) + \ exp[-(^]^(a;,t) 


27r 


27r 


'P 

- expi- 

TT 


Ph'P] + (^DWa^g{x,t) 
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for X G ^Br^,{0) and x = x/r^,. Comparing with (5.5), we find that to leading order 
g, and hence u, satisfy the reactive Robin boundary condition for /3 ^ oo 

DW^g{x,t) ■ X = y^K{^,D) g{x,t) for x G 


where 


K{ip,D) = 


' ^ exp [-(/?] 

Z TT 


1 r erfc(v^) 

1 - 1/ - exp[-^]-^ 


(5.6) 


As —>■ oo, the denominator of K approaches one. This suggests the scaling 




(5.7) 


so that 


as /3 —>■ c». 

We then find g{x^ t) satisfies the desired Robin boundary condition, 

D VxQ{x, t) ■ X = K g{x, t) for x G dB^^^ (0), 

in the limit /3 —)■ oo. More generally, we could impose the equality ^/]3 K(jp, D) = K 
in equation (5.6). Solving for /3, we obtain the following relation between Tp and /3: 

/? = ^ exp[^] - erfc(x/^) exp[^]^ . (5.8) 

In Figure 5.1(a), we compare both formulae (5.7) and (5.8) for AT = £> = 1. As 
expected, they are equivalent in the limit /3 -G oo. In the following section, we will 
present illustrative simulations, confirming that (5.8) provides a more accurate ap¬ 
proximation of the limiting Robin boundary condition. Further improvements to these 
formulas could presumably be made by incorporating a more detailed representation 
of the kinetic boundary layer near i9i3rb(0)- 

5.1. A numerical example illustrating limit /3 —oo. 

As in Section 4.1, we consider the LD model (1.1) for d = 1 in the interval 
n = [0, L] with the linear potential (3.3) near the left-hand boundary only. We 
choose a small time step At and compute the position X{t + At) and the velocity 
V{t + At) from the position X{t) and the velocity V{t) by (4.9)-(4.10). We implement 
adsorbing boundaries at both ends a: = 0 and a; = L of the simulation domain [0, L] 
by terminating the computed trajectory whenever X(t) <0 or X{t) > L. 

In this section, we want to show that both equations (5.7) and (5.8) correctly 
recover the Robin boundary condition in the limit /? —> oo. In particular, we choose 
the values of other parameters equal to 1, namely (compare with (4.11)) 

K = D = L = 1. (5.9) 

We vary /? and we use either equation (5.7) or equation (5.8) to calculate the corre¬ 
sponding value of ^ (these values are plotted in Figure 5.1(a)). 
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Fig. 5.1: (a) Dependence ofTp on (3 computed by equation (5.7) (blue solid line) and 
equation (5.8) (red dashed line) for K = D = 1. 

(b) Splitting probability of escaping through the left boundary as a function (3. In each 
simulation, we estimate the splitting probability as an average over 10® realizations. 
We use £ = 10“®, D = L = K = 1 and tp computed according to equation (5.7) (blue 
squares) or equation (5.8) (red circles). The dashed line denotes the theoretical BD 
result (5.10). 


For each value of (3, we compute 10® trajectories according to (4.9)-(4.10), starting 
from the middle of the domain, i.e. X(0) = F/2, with initial velocity F(0) sampled 
from the normal distribution with zero mean and variance /3D. Each trajectory is 
calculated until it leaves the domain [0, L] either through the left or right boundary 
point. In Figure 5.1(b), we plot the probability that a trajectory leaves the domain 
n through the left boundary (the so-called splitting probability), estimated as the 
fraction of all trajectories which are terminated because X(t) < 0. 

Our goal is to illustrate that equations (5.7) and (5.8) can be used to connect 
the LD model with the limiting Robin boundary problem (2.9)-(2.10). Let n(a:) be 
the probability that the BD particle leaves the domain O = [0, L] through the left 
boundary. Since 


-r^(a:) = 0, for a; e O, 
da;^ 

we find 


with D^(0) =iF(n(0)-l) 
ax 


and n(L) = 0, 


n(x) 


K{L-x) 
D + KL ' 


Since all trajectories start from the middle of the domain, X(0) = Ll2, we have 

for the parameter values given by (5.9). This value is plotted in Figure 5.1(b) as 
the black dashed line. We confirm that the results estimated from simulations ap¬ 
proach (5.10) as /3 —)■ 00 . We also confirm that simulations based on the higher-order 
approximation (5.8) converge more quickly to the limiting Robin boundary problem 
than the simulations based on equation (5.7). 
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6. Discussion. We have considered three parameters, e (potential width), Tp 
(potential height) and (3 (friction constant), and studied several limits of these pa¬ 
rameters which lead to the Robin (reactive) boundary condition (2.10). Parameters 
£ and ^ are shared by both the BD and LD models. In Section 3, we have shown 
that the BD model can recover the Robin boundary condition in the limit £ —0 
and ^ 00 when these parameters are related by (3.13). For the case of the linear 

potential (3.3) this relation can be rewritten as 


= . ( 6 . 1 ) 

The LD model has an additional parameter j3. In Section 5, we have derived two 
formulae (5.7) and (5.8) which relate the LD model to the BD model with the Robin 
boundary condition (2.10). Both results, (5.7) and (5.8), are equivalent in the limit 
P —)■ oo. Equation (5.8) is more accurate for finite values of /?, while (5.7) is simpler 
and easier to interpret. It is given as 


Lp = In 




( 6 . 2 ) 


Considering that (experimentally determinable) parameters K and D are given con¬ 
stants, we can compare our BD result (6.1) with our LD result (6.2). They can be 
both used to specify the height of the potential barrier (^, which is given as a function 
of e in (6.1) and as a function of /3 in (6.2). 

On the face of it, the LD result (6.2) does not depend on the parameter e. How¬ 
ever, the derivation of the LD result (6.2) is only valid for small e. More precisely, 
our conclusion for the LD model (1.1) can be stated as follows: 

(1) If £ <C 1, then Tp is independent of e and can be written in the 

form (6.2); 

(2) If 1 ;:^ £ ^ then Tp is independent of /3 and is given by (6.1). 

In our illustrative computations in Figure 5.1(b) we have used £ = 10“^. In this case, 
the BD result (6.1) implies that Tp = 9.118. We observe that this value is higher 
than the values of Tp plotted in Figure 5.1(a) which are used for our simulations in 
Figure 5.1(b). We can also substitute this value Tp = 9.118 into (6.2). We obtain 
P = 5.224 X 10®. For these values of P and e we are in case (2), so that the LD 
result (6.2) would no longer be applicable. We instead recover the limiting Robin 
boundary condition (in the limit /3 —>■ oo with e = 10“®) by using the BD result (6.1). 

The above results (6.1)-(6.2) can be used to connect the experimentally deter¬ 
minable parameters K and D with parameters of computer simulations to design 
reaction-diffusion models based on LD, i.e. to simulate diffusion-limited bimolecular 
reactions between Kramers particles. There are also other situations where our anal¬ 
ysis will be applicable. One of them is adsorption to surfaces [15,17]. Our set up 
includes interactions with a reactive surface of a sphere and can be used in modeling 
interactions of small molecules with large reactive spheres, for example, for adsorp¬ 
tion of polymers to a surface of a virus [17] or for coating of spherical particles by 
reactive polymers [33] . We also expect our results should be easy to extend to general 
non-spherical reaction surfaces, assuming sufficient regularity. Another possible ap¬ 
plication area is modeling excluded-volume effects. We have observed that the short 
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range repulsive interaction potential (2.6) leads to zero Neumann boundary condi¬ 
tion (3.14) if £ is chosen to approach zero slower than (3.13) (in the limit Tp —>■ oo). A 
similar potential mechanism has been used to enforce Neumann boundary conditions 
on global domain boundaries in [5] , and can be used to model excluded-volume effects 
in models of intracellular macro-molecular crowding [4,6] . 

Since LD requires a smaller time step than the overdamped BD model, the numer¬ 
ical simulations are in general more computationally intensive for LD. However, for 
many biological applications the LD model is only required close to a reactive surface 
(where we have a non-zero potential). In particular, one could replace the compu¬ 
tationally intensive LD model with the BD model in the part of the computational 
domain which is far from the reactive surface [12]. In some applications, one could 
further substitute the BD simulation algorithms by even coarser and more efficient 
simulation techniques, including lattice-based simulations [11,19] or even mean-field 
equations [15,20]. In this way, one could design LD simulation methods which simu¬ 
late intracellular processes on comparable time scales as the BD simulation packages 
which are available in the literature [3,30,34]. In some applications, one could also 
design an efficient first-passage-time scheme by replacing discretization (4.9)-(4.10) 
(which uses fixed time step At) by estimating the next time when a Kramers particle 
becomes sufficiently close to the reactive target [21] . Similar approaches have been 
used to accelerate BD simulations in the literature [27,28,35]. 

Appendix A. A simple example illustrating how specular reflection 
arises in a classical mechanical system. In this section we illustrate with a simple 
example how the specular reflection boundary condition (4.7) arises in a classical 
mechanical system in the absence of noise. 

We consider the simple case of a particle moving in the infinite one-dimensional 
potential. 


kBTp){x) 


-'^x, x<0, 

0, X > 0, 


and experiencing friction, with friction constant p. Note, for the purposes of this 
section we only make use of the value of the potential for x G (—c»,0]. Newton’s 
equation when the particle’s position, X(t), satisfies A(f) < 0 with negative velocity 
is then 


dV 


m—— = 

dt 


—mPV + 


k^Tip 

e 


Using the Einstein Relation this reduces to 


d^ 

dt 


-pv + 


DPip 


the one-dimensional analogue of the LD model (1.1) with the noise term neglected. 
We consider the initial conditions 


U(0) = uo < 0, A(0) = 0, 

and ask for what initial velocity range the molecule successfully moves a distance e 
to the left before changing direction and falling down the potential gradient. This 
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is consistent with a molecule reaching the reactive boundary in the Kramer’s equa¬ 
tion model (2.7). Given the molecule’s negative initial velocity, until the molecule’s 
direction of motion reverses and the molecule moves back across the origin we find 

V{t) = exp [-pt] + 

The velocity then becomes zero at time 

^ In (1 - a), 


where 


a = 


£Vo 

Dp) 


< 0 . 


At this time, the molecule’s position is 

= ^ (« - ln(l - a)) ■ 

The molecule then successfully travels further than —£, analogous to penetrating the 
“reactive boundary” at x = —e, if 


X{t*) < -e. 


or equivalently that 


Dip 

eP 


(a — ln(l 


a)) < —e. 


As e —>• 0 


X{t*) 


-£^o 

2ifPD' 


We therefore find that the molecule successfully penetrates the reactive boundary in 
the limit that e —>■ 0 if and only if 

Wo < -y/2pPD, 


consistent with the specular reflection boundary condition (4.7). 
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